      program visr_n
c     Program to smooth velocity solution using exponential decay distance 
c       function. August, 2024. 
c --------------------------------------------------------------------
      implicit real*8 (a-h,l,o-z)
      parameter (maxsit=150000,maxdat=2500,maxc=50)
      parameter (maxarea=10,maxexcl=20)
      character*100 head
      character*40 gpsfh,gpsfv,outf,crpf,resf,orbf
      character*40 gps_h_outf,gps_v_outf,gps_prdf,area_excf
      character*8 stnlh(maxsit+1)
      integer*8 nstn,nstnh,nstnv,ncrp,nsite,ms,nold,ndata
      dimension area_lon(maxarea,maxexcl),area_lat(maxarea,maxexcl)
      dimension ar_lon(maxexcl),ar_lat(maxexcl)
      dimension i_area(maxexcl)
      dimension lonlh(maxsit+1),latlh(maxsit+1)
      dimension lonh(maxsit),lath(maxsit),uxlh(maxsit),sxlh(maxsit),
     .  uylh(maxsit),sylh(maxsit)
      dimension glon(maxsit),glat(maxsit),qlon(1),qlat(1)
      dimension posh(2,maxsit)
      dimension areah(maxsit)
      dimension careah(maxsit)
      dimension alon(maxc),alat(maxc),blon(maxc),blat(maxc)
      dimension bbz(6),dcs(maxc),dsn(maxc),tx(1),ty(1)
      dimension a(maxdat,maxdat),aat(maxdat,maxdat),indx1(maxdat)
      dimension bb(6),ainv(6,6),b(maxdat),bbt(maxdat)
      dimension prd0(90),asv(90,3),bsv(90)
      dimension vxs(maxsit+1),vys(maxsit+1),vx_sig(maxsit+1),
     1    vy_sig(maxsit+1),vzs(maxsit+1),vz_sig(maxsit+1),tmp(maxsit+1)
      integer*8 m,idat
      common /veldath/nstnh,lonh,lath,uxlh,sxlh,uylh,sylh,uzlh,szlh,
     1     cxyh,areah
      common /crpdat/ncrp,alon,alat,blon,blat,dcs,dsn
      common /ainv_vh/ainv,bb,rtau,wtt,chisq,nslct
      common /prmtr/pi,cov,cutoff_dis,wt_az
*
*   setup the constants                          
*
      pi = 4.0d0*datan(1.0d0)
      cov = pi/180.0d0
      ms = 1
      is1 = 0
*
*   read input file names and options 
*
      read(5,'(a)') gpsfh                                      ! input horizontal gps data file
      read(5,'(a)') outf                                       ! output file
      read(5,*) rtau_h                                         ! default smoothing distances 
      read(5,*) is_wght                                        ! distance weighting scheme: 1=gaussian, 2=quadratic
      read(5,*) id_wght                                        ! spatial weighting scheme: 1=azimuth, 2=voronoi area

      open(4,file=gpsfh,status='old')                           ! open horizontal gps velocity data
      open(41,file=outf,status='unknown')                       ! open output file
      write(41,24)
24    format("  Longitude  Latitude      Ve       Vn      See       Snn
     1      Sen       W")
      read(4,*) 
      i = 1
30    read(4,35,end=40) lonlh(i),latlh(i),uxlh(i),sxlh(i),
     1   uylh(i),sylh(i)
      if (lonlh(i).gt.180.0d0 ) lonlh(i) = lonlh(i)-360.0d0
c      if (sxlh(i).lt.rsgah) sxlh(i)=rsgah  
c      if (sylh(i).lt.rsgah) sylh(i)=rsgah  
      i = i + 1
      goto 30
35    format(2f9.3,4f9.2)
40    nstnh=i-1
      close(4)
c
      if (nstnh.gt.maxsit) then            ! check limitation of data array input
        write(6,62) maxsit
        stop
      endif
62    format('# of stations in input file > maxsit',i5,
     .   ', please increase maxsit ...')
*
*   compute mass center of network
*
      xmean=0.0d0
      ymean=0.0d0
      do i=1,nstnh
        xmean=xmean+lonlh(i)
        ymean=ymean+latlh(i)
      enddo
      xmean=xmean/nstnh
      ymean=ymean/nstnh
*
*  convert geodetic coordinates into cartesian coordinates
*
      do i=1,nstnh
        lonh(i)=lonlh(i)
        lath(i)=latlh(i)
      enddo
      call llxy(ymean,xmean,lath,lonh,nstnh)
c
      read(5,*) ncrp
      if (ncrp.gt.maxc) then
        write(6,61) maxc
        stop
      endif
61    format('# of creep faults > maxc',i3,
     .   ', please increase maxc ...')
      if (ncrp .ne. 0) then
        read(5,'(a)') crpf
        open(4,file=crpf,status='old')                   ! open file for creep fault parameter input
        do i = 1, ncrp
          read(4,*) alon(i),alat(i),blon(i),blat(i)
          if (alon(i).gt.180.0d0 ) alon(i) = alon(i)-360.0d0
          if (blon(i).gt.180.0d0 ) blon(i) = blon(i)-360.0d0
        enddo
      endif
*
*  convert geodetic coordinates into cartesian coordinates for creep fault and
*  compute creep fault's azimuth
*
      if ( ncrp.ne.0 ) then 
        call llxy(ymean,xmean,alat,alon,ncrp)
        call llxy(ymean,xmean,blat,blon,ncrp)
        do i = 1, ncrp
          dx = blon(i) - alon(i)
          dy = blat(i) - alat(i)
          ds = dsqrt(dx**2 + dy**2)
          dcs(i) = dx/ds
          dsn(i) = dy/ds
        enddo
      endif
*
*   setup constants for data reweighting according to options                          
*
      if (is_wght .eq. 1) cutoff_dis = 5.15d0
      if (is_wght .eq. 2) cutoff_dis = 3.0*rtau_h
      wt_az = 0.25d0                                                         ! relative weight of mean azimuth to individual azimuth
c
c     gps site horizontal velocity interpolation
c
      do i = 1, nstnh
         yy = lath(i)
         xx = lonh(i)
         call visr_h(yy,xx,rtau_h,
     1      is_wght,id_wght,indxx,iout,is1)                   !  horizontal gps velocity interpolation and uncerstainty estimate
         write(41,28) lonlh(i),latlh(i),bb(1),bb(2),bb(3),
     1    bb(5),bb(4),bb(6)                                   !  output horizontal GPS velocity result from bootstrapping
      enddo
28    format(2f10.3,2f9.3,4f10.3)
c
      stop
      end
c
c     program to judge whether a site is inside a enclosed region.  9/4/2021.
c
      subroutine inside(nd,np,xp,yp,xsite,ysite,insd)
      implicit real*8 (a-h, o-z)
      real*8 xp(np),yp(np),azi(np)
      cov = 3.14158/180.0

      do i = 1, nd
        dx = xp(i) - xsite
        dy = yp(i) - ysite
        azi(i) = datan2(dy,dx)/cov
      enddo
      azi_sum = 0
      do i = 1, nd-1
         dazi = azi(i+1) - azi(i)
         if (dazi .gt. 180.0) dazi = dazi - 360.0
         if (dazi .lt. -180.0) dazi = dazi + 360.0
         azi_sum = azi_sum + dazi
      enddo
      if (abs(azi_sum) .gt. 180.0) then
         insd = 1
      else
         insd = 0
      endif
      return
      end
c
      subroutine llxy(slatm,slonm,slat,slon,m)
*
*   subroutine to transfer latitude and longitude into local x y.
*
*   input:
*      slatm: latitude of reference point for local coordinates, in degree.
*      slonm: longitude of reference point for local coordinates, in degree.
*      slat : array coordinate of latitude to be transferred, in degree.
*      slon : array coordinate of longitude to be transferred, in degree.
*      m    : number of points to be transferred.
*   output:
*      slat : y array after transformation in km.
*      slon : x array after transformation in km.
*
      implicit real*8 (a-h,o-z)
      integer*8 m
      dimension t(3,3),vp(3),v(3)
      dimension slat(m),slon(m)
*
      pi = 4.0d0*datan(1.0d0)
      rlatc = slatm*pi/180.0d0
      rlonc = slonm*pi/180.0d0
*
*     calculate local radius of curvature (r) using reference earth nad 
*     1983 (semi-major axis is 6378.137,flattening factor is 1/298.2572)
*
      flat = 1.0d0/298.2572d0
      esq = 2.0d0*flat - flat**2
      q = 1.0d0 - esq*dsin(rlatc)*dsin(rlatc)
      r = 6378.137d0*dsqrt(1.0d0 - esq)/q
*
*     construct transformation matrix
*
      t(1,1) =  dsin(rlatc)*dcos(rlonc)
      t(1,2) =  dsin(rlatc)*dsin(rlonc)
      t(1,3) = -dcos(rlatc)
      t(2,1) = -dsin(rlonc)
      t(2,2) =  dcos(rlonc)
      t(2,3) =  0.0d0
      t(3,1) =  dcos(rlatc)*dcos(rlonc)
      t(3,2) =  dcos(rlatc)*dsin(rlonc)
      t(3,3) =  dsin(rlatc)
*
*     calculate xl,yl,dis (vector in local coordinate)
*
      do 9 i = 1,m
        slat(i) = slat(i)*pi/180.0d0
        slon(i) = slon(i)*pi/180.0d0
        v(1) = dcos(slat(i))*dcos(slon(i))
        v(2) = dcos(slat(i))*dsin(slon(i))
        v(3) = dsin(slat(i))
        do j = 1,2
          vp(j) = 0.0d0
          do k = 1,3
            vp(j) = vp(j) + t(j,k)*v(k)
          end do
        end do
        slon(i) =  r*vp(2)
        slat(i) = -r*vp(1)
9     continue

      return
      end

*  Subroutine to conduct LU back substitution
*
      subroutine dlubksb(a,n,np,indx,b)
      implicit real*8 (a-h,o-z)
      dimension a(np,np),indx(np),b(np)
      ii=0
      do 12 i=1,n
        ll=indx(i)
        sum=b(ll)
        b(ll)=b(i)
        if (ii.ne.0)then
          do 11 j=ii,i-1
            sum=sum-a(i,j)*b(j)
11        continue
        else if (sum.ne.0.0d0) then
          ii=i
        endif
        b(i)=sum
12    continue
      do 14 i=n,1,-1
        sum=b(i)
        if(i.lt.n)then
          do 13 j=i+1,n
            sum=sum-a(i,j)*b(j)
13        continue
        endif
        b(i)=sum/a(i,i)
14    continue
      return
      end
 
*  Subroutine to conduct LU decomposition
*
      subroutine dludcmp(a,n,np,indx,d)
      implicit real*8 (a-h,o-z)
      parameter (nmax=50000,tiny=1.0e-20)
      dimension a(np,np),indx(np),vv(nmax)
      d=1.0d0
      do 12 i=1,n
        aamax=0.0d0
        do 11 j=1,n
          if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j))
11      continue
        if (aamax.eq.0.0d0) write (6,100) i
        vv(i)=1.0d0/aamax
12    continue
100   format('singular matrix, i= ',i4)

      do 19 j=1,n
        if (j.gt.1) then
          do 14 i=1,j-1
            sum=a(i,j)
            if (i.gt.1)then
              do 13 k=1,i-1
                sum=sum-a(i,k)*a(k,j)
13            continue
              a(i,j)=sum
            endif
14        continue
        endif
        aamax=0.0d0
        do 16 i=j,n
          sum=a(i,j)
          if (j.gt.1)then
            do 15 k=1,j-1
              sum=sum-a(i,k)*a(k,j)
15          continue
            a(i,j)=sum
          endif
          dum=vv(i)*abs(sum)
          if (dum.ge.aamax) then
            imax=i
            aamax=dum
          endif
16      continue
        if (j.ne.imax)then
          do 17 k=1,n
            dum=a(imax,k)
            a(imax,k)=a(j,k)
            a(j,k)=dum
17        continue
          d=-d
          vv(imax)=vv(j)
        endif
        indx(j)=imax
        if(j.ne.n)then
          if(a(j,j).eq.0.0d0)a(j,j)=tiny
          dum=1./a(j,j)
          do 18 i=j+1,n
            a(i,j)=a(i,j)*dum
18        continue
        endif
19    continue
      if(a(n,n).eq.0.0d0)a(n,n)=tiny
      return
      end
*
*  Subroutine for GPS horizontal velocity interpolation
*
      subroutine visr_h(yy,xx,rtau_h,
     .   is_wght,id_wght,indxx,iout,is1)
      implicit real*8 (a-h,l,o-z)
      integer*8 nstn,ncrp
      parameter (maxsit=150000)
      dimension lon(maxsit),lat(maxsit),uxl(maxsit),sxl(maxsit),
     .          uyl(maxsit),syl(maxsit),cxy(maxsit)
      dimension alon(50),alat(50),blon(50),blat(50),
     .          dcs(50),dsn(50)
      dimension crp_strk1(50),crp_strk2(50),is0(50),icrp(50)
      dimension istn(maxsit),azi(maxsit),wght(maxsit)
      dimension a(2*maxsit,6),b(2*maxsit),aa(6,6),aasv(6,6),
     .          ainv(6,6),indx(6),bb(6),v(6),az(maxsit,3)
      dimension a_s(2*maxsit,6),aa_s(6,6)
      common /veldath/nstn,lon,lat,uxl,sxl,uyl,syl,uzl,szl,cxy,area
      common /crpdat/ncrp,alon,alat,blon,blat,dcs,dsn
      common /ainv_vh/ainv,bb,rtau,wtt,chisq,nslct
      common /prmtr/pi,cov,cutoff_dis,wt_az

      crp_dis = 30               ! criterion distance for exemption of azimuthal coverage for sites near creep fault
      iout = 0
      i_ctc = 0
*
*  compute relationship between interpolation point and creep fault
*
      do i = 1, ncrp

        icrp(i) = 0

        dx = alon(i) - xx
        dy = alat(i) - yy                                ! compute vector V_a1 = F_a - R_int
        dra = sqrt(dx*dx + dy*dy)
        ang1 = datan2(dy,dx)/cov                         ! compute azimuth angle of vector from fault end point a to interpolation spot
        y1 = -1.0d0*(dy*dcs(i)-dx*dsn(i))                ! compute V_1 = V_a1 X F
        dx = blon(i) - xx
        dy = blat(i) - yy                                ! compute vector Vb = Fb - R
        drb = sqrt(dx*dx + dy*dy)
        ang2 = datan2(dy,dx)/cov                         ! compute azimuth angle of fault end point b measured from interpolation point
        if (dra .lt. crp_dis .or. drb .lt. crp_dis) i_ctc = 1
        if ( ang1.eq.ang2 ) then
          is0(i) = 0                                     !  interpolation point is located on the extension of creep fault ?
          goto 20
        endif 
        if ( y1.eq.0.0d0 ) then
          indxx = 4                                      !  interpolation point is located on creep fault
          goto 60  
        endif

        is0(i) = y1/abs(y1)

20      if (ang2 .gt. ang1) then                          ! set angle 1 > angle 2, mesured counter-clockwise
          crp_strk1(i) = ang1
          crp_strk2(i) = ang2
        else
          crp_strk1(i) = ang2
          crp_strk2(i) = ang1
        endif

        dang = crp_strk2(i) - crp_strk1(i)                ! compute angle difference
        if ( dang.gt.180.0d0 ) then
          ang_tmp = crp_strk1(i) + 360.0d0
          crp_strk1(i) = crp_strk2(i)
          crp_strk2(i) = ang_tmp
          icrp(i) = 1
        endif

      enddo
*
*  start data selection.
*
        indxx = 0
        nslct = 0
        rtau = rtau_h

        do 10 i = 1, nstn                             ! iteration to select gps data sites

          dx = lon(i) - xx
          dy = lat(i) - yy
          dr = dsqrt(dy**2+dx**2)
          ang_s1 = datan2(dy,dx)/cov                  ! compute station azimuth viewing from interpolation site
          rt = dr/rtau

          if (rt .le. cutoff_dis) then
            do j = 1, ncrp
              if (icrp(j).eq.1.and.ang_s1.lt.0.d0) ang_s1=ang_s1+360.d0
 
              dx = alon(j) - lon(i)
              dy = alat(j) - lat(i)                  ! compute V_a2 = F_a - Rgps
              y1 = -1.0d0*(dy*dcs(j)-dx*dsn(j))      ! compute V_2 = V_a2 X F 
              
              if ( y1.eq.0.0d0) then                 ! site on fault 
                 is1 = 0
              else
                 is1 = y1/abs(y1)
              endif
             
              if (is1.ne.is0(j)) then             ! gps site and interpoltion spot are on oposite sides of fault
              if (ang_s1.ge.crp_strk1(j).and.
     .          ang_s1.le.crp_strk2(j)) then      ! gps site in shaddow, discard
                goto 10
              endif
              endif

            enddo
            nslct = nslct + 1
            istn(nslct) = i
          endif
10      continue

      if (nslct.lt.3) then
        indxx = 1
        goto 30
      endif                                                                  ! no interpolation for n_slct < 3

      if (nslct.gt.maxsit) then
        write(6,62) maxsit
        stop
      endif

      if (i_ctc .eq. 1) goto 21
62    format('# of selected stations > maxsit',i5,
     .   ', please increase maxsit ( visr_core )...')
*
*  estimate azimuthal coverage.
*
      do i = 1, nslct
        dx = lon(istn(i)) - xx
        dy = lat(istn(i)) - yy
        azi(i) = datan2(dy,dx)/cov
      enddo
      azimax = azi(1)
      azimin = azi(1)
      do i = 2, nslct
        azimax = max(azimax,azi(i))
        azimin = min(azimin,azi(i))
      enddo
      dazim1 = azimax - azimin
      if (dazim1.le.180.0d0) then
        indxx = 2
        goto 30
      endif                                                                  ! no interpolation for dazi < pi

      do i = 1, nslct
        if (azi(i).lt.0.d0 ) azi(i) = azi(i) + 360.0d0
        if ( i.eq.1 ) then
          azimax = azi(i)
          azimin = azi(i)
        else
          azimax = max(azimax,azi(i))
          azimin = min(azimin,azi(i))
        endif
      enddo
      dazim2 = azimax - azimin
      if (dazim2.le.180.0d0) then
        indxx = 2
        goto 30
      endif                                                                  ! no interpolation for dazi < pi

21    iout = 1
*
      npr = 2*nslct
      wtt = 0.0d0
      do i = 1, nslct
        dy = lat(istn(i))-yy
        dx = lon(istn(i))-xx
        dr = dsqrt(dx**2+dy**2)
        if (is_wght .eq. 1) wti = dexp(-dr/rtau)
        if (is_wght .eq. 2) wti = 1/((dr/rtau)**2+1)
        wtt = wtt + 2*wti

        a(2*i-1,1)=wti/sxl(istn(i))
        a(2*i-1,2)=0.0d0
        a(2*i-1,3)=wti*dx/sxl(istn(i))
        a(2*i-1,4)=wti*dy/sxl(istn(i))
        a(2*i-1,5)=0.0d0
        a(2*i-1,6)=wti*dy/sxl(istn(i))
        a(2*i,1)=0.0d0
        a(2*i,2)=wti/syl(istn(i))
        a(2*i,3)=0.0d0
        a(2*i,4)=wti*dx/syl(istn(i))
        a(2*i,5)=wti*dy/syl(istn(i))
        a(2*i,6)=-1.0*wti*dx/syl(istn(i))

        b(2*i-1)=wti*uxl(istn(i))/sxl(istn(i))
        b(2*i)=wti*uyl(istn(i))/syl(istn(i))

      enddo

      if (wtt .lt. wt0) then
        indxx = 3
        goto 30
      endif                                                                  ! no interpolation for total weight < weighting threshold

      ww = 0.0d0 
      do i = 1, npr
         ww = ww + b(i)**2
      enddo

      do i=1,6
         do j=1,i
            aa(i,j)=0.0d0
            do k=1,npr
               aa(i,j)=aa(i,j)+a(k,i)*a(k,j)
            enddo
            aa(j,i)=aa(i,j)
         enddo
      enddo
         
      do i = 1, 6
         do j = 1, 6
            aasv(i,j) = aa(i,j)
         enddo
      enddo

      do i=1,6
         bb(i)=0.0d0
         do j=1,npr
            bb(i)=bb(i)+a(j,i)*b(j)
         enddo
      enddo
*
*  solve b=a*x.
*
      call dludcmp(aa,6,6,indx,d)
      call dlubksb(aa,6,6,indx,bb)                                  ! solve for solution
*
*  invert n.
*
      do i=1,6
         do j=1,6
            ainv(i,j)=0.0d0
         enddo
         ainv(i,i)=1.0d0
      enddo

c      call mat_aba(bb,aasv,6,xaax)
c      chisq = ww - xaax

      goto 60

30    continue

60    return
      end
